arXiv: 1509.00677vl [astro-ph.IM] 2 Sep 2015 


To be submitted to PASP 


Performance analysis of the Least-Squares estimator in Astrometry 

Rodrigo A. Lobos, Jorge F. Silva, 

Departamento de Ingenieria Electrica, Facultad de Ciencias Fisicas y Matemdticas, 
Universidad de Chile, Beauchef 850, Santiago, Chile 

rlobos, josilva Oing.uchile.cl 
Rene A. Mendez, 

Departamento de Astronomia, Facultad de Ciencias Fisicas y Matemdticas, 
Universidad de Chile, Casilla 36-D, Santiago, Chile 

rmendezSu.uchile.cl 
and 

Marcos Orchard 

Departamento de Ingenieria Electrica, Facultad de Ciencias Fisicas y Matemdticas, 
Universidad de Chile, Beauchef 850, Santiago, Chile 

morchardSing.uchile.cl 

ABSTRACT 

We characterize the performance of the widely-used least-squares estimator in as¬ 
trometry in terms of a comparison with the Cramer-Rao lower variance bound. In this 
inference context the performance of the least-squares estimator does not offer a closed- 
form expression, but a new result is presented (Theorem]^ where both the bias and the 
mean-square-error of the least-squares estimator are bounded and approximated ana¬ 
lytically, in the latter case in terms of a nominal value and an interval around it. From 
the predicted nominal value we analyze how efficient is the least-squares estimator in 
comparison with the minimum variance Cramer-Rao bound. Based on our results, we 
show that, for the high signal-to-noise ratio regime, the performance of the least-squares 
estimator is significantly poorer than the Cramer-Rao bound, and we characterize this 
gap analytically. On the positive side, we show that for the challenging low signal- 
to-noise regime (attributed to either a weak astronomical signal or a noise-dominated 
condition) the least-squares estimator is near optimal, as its performance asymptoti¬ 
cally approaches the Cramer-Rao bound. However, we also demonstrate that, in general. 
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there is no unbiased estimator for the astrometric position that can precisely reach the 
Cramer-Rao bound. We validate our theoretical analysis through simulated digital- 
detector observations under typical observing conditions. We show that the nominal 
value for the mean-square-error of the least-squares estimator (obtained from our the¬ 
orem) can be used as a benchmark indicator of the expected statistical performance of 
the least-squares method under a wide range of conditions. Our results are valid for an 
idealized linear (one-dimensional) array detector where intra-pixel response changes are 
neglected, and where flat-fielding is achieved with very high accuracy. 

Subject headings: Astrometry, parameter estimation, Cramer-Rao bound, Least-Squares 
estimator, performance analysis. 


Motivation 


Astrometry, that branch of observational astronomy that deals with the precise and accurate 
estimation of angular positions of light-emitting (usually point-like) sources projected against the 


celestial sphere, is the oldest technique employed in the study of the heavens (Hpg 2009 Reffert 


2009; H0g||2OlT ). Repeated measurements of positions, spread over time, allow a determination of 
the motions and distances of theses sources, with astrophysical implications on dynamical studies 
of stellar systems and the Milky Way as a whole. With the advent of solid-state detectors and all- 
digital techniques applied to radio-interferometry and specialized ground- and space-based missions, 
astrometry has been revolutionized in recent years, as we have entered a high-precision era in 
which this technique has started to play an increasingly important role in all areas of astronomy, 
astrophysics (van Altena 2013), and cosmology ( Lattan^|2012 ). 

Current technology, based on two-dimensional discrete digital detectors (such as charged cou¬ 
pled devices - CCDs), record a (noisy) image (on an array of photo-sensitive pixels) of celestial 
sources, from which it is possible to estimate both their astrometry and photometry, simultane¬ 


ously (Howell 2006). The inference problem associated to the determination of these quantities is 


at the core of the astrometric endeavor described previously. 

A number of techniques have been proposed to estimate the location and flux of celestial 
sources as recorded on digital detectors. In this context, estimators based on the use of a least- 
squares error function (LS hereafter) have been widely adopted ( Stetson||1987 King||1983 Alard &: 
Lupton||1998 Honeycutt||1992 Cameron et al.||2006 ). The use of this type of decision rule has been 
traditionally justified through heuristic reasons and because they are conceptually straightforward 
to formulate based on the observation model of these problems. Indeed, the LS approach was 


the classical method used when the observations were obtained with analog devices (van Altena 


&: Auer 1975; Auer & van Altena 1978) (which corresponds to a Gaussian noise model for the 
observations, different from that of modern digital detectors, which is characterized instead by a 
Poisson statistics) and, consequently, the LS method was naturally adopted from the analogous to 
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the digital observational setting 

In contemporary astrometry (Gaia, for instance), stellar positions will be obtained by optimiz¬ 
ing a likelihood function (see, e.g., |Lindegren| ( [2000 ), which uses the equivalent of our equations 
and (13) in Sections |2 . 1| and [3(T] respectively), not by LS. Nevertheless, since LS methods offer com¬ 


putationally efficient implementations and have shown reasonable performance (Lee & van Altena 


1983; Stone||1989) they are still widely used in astrometry either on general-purpose software pack¬ 


ages for the analysis of digital images such as DAOPHOT ( [Stetson 1987), or on dedicated pipelines, 
such as that adopted in the Sloan Digital Sky Survey survey (SDSS hereafter, Lupton et ah (2001); 


Lupton (2007)). For example, in DAOPHOT astrometry (and photometry) are obtained through 


a two-step process which involves a LS minimization of a trial function (e.g., a bi-dimensional 
Gaussian, see Stetson (1987 equation 6), equivalent to our one-dimensional case in equation (16), 
Section 3.1), and then applying a correction by means of an empirically determined look-up table 
(also computed performing a LS on a set of high signal-to-noise ratio images distributed over the 
field of view of the image, see Stetson (1987, equation 8)). This last step accounts for the fact that 


the PSF of the image under analysis may not be exactly Gaussian. The SDSS pipeline (Lupton 


et al. (2001)) obtains its centroids also through a two-step process: First it fits a Karhunen-Loeve 


m 


transform (KL transform hereafteij^ see, e.g., Najim (2008)) to a set of isolated bright stars i 
the field-of-view of the image, and then it uses the base functions determined in this way, to fit the 
astrometry and photometry for the object(s) under consideration using a LS minimization scheme 
(see Lupton (2007, equation (5))). Both codes, DAOPHOT and the SDSS pipeline have been ex¬ 
tensively used and tested by the astronomical community, giving very reliable results (see, e.g., 


Pier et al. (2003); Becker et al. (2007)). 


Considering that LS methods are still in use in astrometry, and driven by the increase in the 
intrinsic precision available by the new detectors and instrumental settings, and by the fact that 
CCDs will likely continue to be the detector of choice for the focal-plane in science-quality imaging 
application at optical wavelengths for both space- as well as ground-borne program^ it is timely to 
re-visit the pertinence of the use of LS estimators. Indeed, in the digital setting, where we observe 
discrete samples (or counts) on a photon integrating device, there is no formal justification that the 
LS approach is optimal in the sense of minimizing the mean-square-error (MSE) of the parameter 
estimation, in particular for astrometry, which is the focus of this work. 


The question of optimality (in some statistical sense) has always been in the interest of the 
astronomical community, in particular the idea of characterizing fundamental performance bounds 
that can be used to analyze the efficiency of the adopted estimation schemes. In this context, 


^See footnote 11 in Section 


0 


^Also referred to as “principal component analysis”, “proper orthogonal decomposition”, “empirical orthogonal 
functions” (a term used in meteorology and geophysics), or the Hotelling transform, see Gray & Davisson ( 2010| . 

^See, e.g., http://www.bnl.gov/cosmo2013/ 
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we can mention some seminal works on the nse of the celebrated Cramer-Rao (CR hereafter) 
bound in astronomy by Lindegren (1978); Jakobsen et al. ( 1992[); [Zaccheo et al.| (1995); Adorf 


(1996); and Bastian (2004). The CR bound is a minimum variance (MV) bound for the family of 


unbiased estimators (Rao 1945 Cramer 1946). In astrometry and photometry this bound has offered 


meaningful closed-form expressions that can be used to analyze the complexity of the inference task, 
and its dependency on key observational and design parameters such as the position of the object in 
the array, the intensity of the object, the signal-to-noise ratio (SNR hereafter), and the resolution of 
the instrument ( Winick||1986 Perryman et al.||1989 Lindegren||2010 Mendez et al.||2013 2014). In 


particular, for photometry, Perryman et al. (1989) used the CR bound to show that the LS estimator 


is a good estimator, achieving a performance close to the limit in a wide range of observational 
regimes, and approaching very closely the bound at low SNR. In astrometry, on the other hand. 


Mendez et al. (2013, 2014) have recently studied the structure of this bound and have analyzed 


its dependency with respect to important observational parameters, under realistic astronomical 
observing conditions. In those works, closed-form expressions for the bound were derived in a 
number of important settings (high spatial resolution, low and high SNR), and their trends were 
explored across angular resolution and the position of the object in the array. As an interesting 
outcome, the analysis of the CR bound allows us to find the optimal pixel resolution of the array 
for a given setting, as well as providing formal justification to some heuristic techniques commonly 


used to improve performance in astrometry, like dithering for undersampled images (Mendez et al 


2013 Sec. 3.3). 


The specific problem of evaluating the existence of an estimator that achieves the CR bound 


has not been covered in the literature, and remains an interesting open problem. On this, Mendez 


et al. (2013) have empirically assessed (using numerical simulations) the performance of two LS 


methods and the maximum-likelihood (ML hereafter) estimator, showing that their variances fol¬ 
low very closely the CR limit in some specific regimes. In this paper, we analyze in detail the 
performance of the LS estimator with respect to the CR bound, with the goal of finding concrete 
regimes, if any, where this estimator approaches the CR bound and, consequently, where it can be 
considered an efficient solution to the astrometric problem. This application is a challenging one, 
because estimators based on a LS type of objective function do not have a closed-form expression 
in astrometry. In fact, this estimation approach corresponds to a non linear regression problem, 
where the resulting estimator is implicitly defined. As a result, no expressions for the performance 
of the LS estimator can be obtained analytically. To address this issue, our main result (Theorem 
Section 3.2) derives expressions that bound and approximate the variance of the LS estimator. 


Our approach is based on the work by So et al. (2013), where the authors tackle the problem of 


approximating the bias and MSE of general estimators that are the solution of an optimization 


problem. In Fessler (1996) another methodology is given to approximate the variance and mean 


of implicitly defined estimators, which has been applied to medical imaging and acoustic source 


localization (Raykar et al. 2005). 


The main result of our paper is a refined version of the result presented in So et al. (2013), where 
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one of their key assumptions, which is not applicable in our estimation problem, is reformulated. 
In this process, we derive lower and upper bounds for the MSE performance of the LS estimator. 
Using these bounds, we analyze how closely the performance of the LS estimator approaches the CR 
bound across different observational regimes. We show that for high SNR there is a considerable 
gap between the CR bound and the performance of the LS estimator. Remarkably, we show that for 
the more challenging low SNR observational regime (weak astronomical sources), the LS estimator 
is near optimal, as its performance is arbitrarily close to the CR bound. 

The paper is organized as follows. Section introduces the problem, notation, as well as some 
preliminary results. Section represents the main contribution, where Theorem and its inter¬ 
pretation are introduced. Section shows numerical analyses of the performance of LS estimator 
under different observational regimes. Finally Section provides a summary of our results, and 
some final remarks. 


2. Introduction and preliminaries 

In this section we introduce the problem of astrometry as well as concepts and definitions that 
will be used throughout the paper. For simplicity, we focus on the 1-D scenario of a linear array 
detector, as it captures the key conceptual elements of the problenij^ 


2.1. Astrometry and photometry based on a photon integrating device 

The specific problem of interest is the inference of the position of a point source. This source 
is parameterized by two scalar quantities, the position of the object Xc £ R in the arrajQ and 
its intensity (or brightness, or flux) that we denote by F G M"*". These two parameters induce a 
probability p over an observation space that we denote by X. More precisely, given a point 
source represented by the pair (xc, F), it creates a nominal intensity profile in a photon integrating 
device, typically a CCD, which can be generally written as: 


= F -Hx- Xc,a), 


( 1 ) 


where 4>{x — Xc, cr) denotes the one dimensional normalized point spread function (PSF) evalu¬ 
ated on the pixel coordinate x — Xc, and where fi is a generic parameter that determines the width 
(or spread) of the light distribution on the detector (typically a function of wavelength and the 
quality of the observing site, see Section]^ (see Mendez et al. (2013, 2014) for more details). 


"^This analysis can be extended to the 2-D case as shown in 

®This is of course related to an angular position in the sky, measured in arcsec, through the “plate-scale”, which 
is an optical design feature determined by the instrument plus telescope configuration. 


Mendez et al. 


(2013 
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The profile in equation Q is not observed directly, but through three sources of perturba¬ 
tions: First, an additive background which accounts for the photon emissions of the open (dif¬ 
fuse) sky and the contributions from the noise of the instrument itself (the read-out noise and 

modeled by Bi in equation Second, an intrinsic 


dark-current (Gilliland 


1992 


Tyson 


1986 


uncertainty between the aggregated intensity (the nominal object brightness plus the background) 
and the actual measurements, denoted by li in what follows, which is modeled by independent 
random variables that obey a Poisson probability law. And, finally, we need to consider the spatial 
quantization process associated with the pixel-resolution of the detector as specified by gi{xc) in 
equations Q and ([sjj^ Including these three effects, we have a countable collection of indepen¬ 
dent and not identically distributed random variables (observations or counts) {/j : i G Z}, where 
li ~ Poisson{\i{xc, F)), driven by the expected intensity at each pixel element i, given by: 


Xi{xc, F) = E{/j} = F ■ gi(xc) +Bi, Vi G Z 


( 2 ) 


=F{x,,F) 


and. 


rXi+A.x/2 


gi{xc)^ / (l){x - Xc,a) dx, yi e Z, (3) 

J Xi — Axl2 

where E represents the expectation value of the argument, and {xi : i G Z} denotes the standard 
uniform quantization of the real line-array with resolution Ax > 0, i.e., Xj+i — xt = Ax for all 
i G Z. In practice, the detector has a finite collection of measured elements (or pixels) xi, ...,Xn, 
then a basic assumption here is that we have a good coverage of the object of interest, in the sense 
that for a given position Xc- 


n 

'^gi{xc) ~ '^gi{xc) = / 

i=l iGIa ~ 


4>{x — Xc, cr) dx = 1. 


( 4 ) 


Note that equation Q adopts the idealized situation where every pixel has the exact same response 
function (equal to unity), or, equivalently, that our flat-field process has been achieved with minimal 
uncertainty. It also assumes that the intra-pixel response is uniform. The latter is more important 


in the severely undersampled regime (see, e.g., Adorf (1996, Figure 1)) which is not explored in this 


paper. However a relevant aspect of data calibration is achieving a proper flat-helding which can 
affect the correctness of our analysis and the form of the adopted likelihood function (see below). 

At the end, the likelihood of the joint random observation vector = (/i,.., In) (with values 
in N”"), given the source parameters {xc,F), is given by: 


L(r; x„ F) = ^)(/2) • • • fxc.ixcF){In), Vr G 


( 5 ) 


®Note that pixel-convolved point-spread function gi{xc) is sometimes referred to as the the ’’pixel response func¬ 
tion” . 
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where fx{x) = 




denotes the probability mass function (pmf) of the Poisson law (Gray & 


Davisson [2010 )■ We emphasize that equation ([^ assumes that the observations (contained in the 


individual pixels, denoted by the index i) are independent (although not identically distributed, 
since they follow A*). Of course, this is only an approximation to the real situation since it implies, 
in particular, that we are neglecting any electronic defects or features in the device such as, e.g., the 


cross-talk present in multi-port CCDs ( 

Freyhammer et al. 

2001 

), or read-out correlations such as 

the odd-even column effect in IR detectors (Mason ||2007 

), as well as calibration or data reduction 

deficiencies (e.g., due to inadequate flat-fielding ( 

Gawiser et al. 112006)) that may alter this idealized 


detection process. A serious attempt is done by manufacturers and observatories to minimize the 
impact of these defects, either by an appropriate electronic design, or by adjusting the detector 
operational regimes (e.g., cross-talk can be reduced to less than 1 part in 10^ by adjusting the 
readout speed and by a proper reduction process, Freyhammer et al. ( 2001[ )). In essence, we 
are considering an ideal detector that would satisfy the proposed likelihood function given by 
equation ([^, in real detectors the likelihood function could be considerably more complej^^ 


We can formulate the astrometric and photometric estimation task, as the problem of char¬ 
acterizing a decision rule r^O : N”' —>• 0, with 0 = being a parameter space, where given an 
observation I” the parameters to be estimated are In other words, Tn{I^) 

gives us a prescription (or statistics) that would allow us to estimate the underlying parameters 
{xc, F) (the estimated parameters are denoted by {xc, F)) from the available data vector I”. 


In the simplest scenario, in which one is interested in determining a single (unknown) parameter 
9 (e.g., in our case either Xc or F, assuming that all other parameters are perfectly well known), a 
commonly used decision rule adopted in statistics to estimate this parameter (the estimation being 
6) is to consider the prescription of minimum variance (denoted by Var), given by: 


9{) = arg min I/ar (a (/"■)) 

a 

= argminE7-n...,/n («(-^’^) “ ^)^ (6) 

where “argmin” represents the argument that minimizes the expression, while a is a generic 
variable representing the parameter to be determined. Note that in the last equality we have 
assumed that 9{) is an unbiased estimator of the parameter (i.e., that E(0) = 9), so that under 
this rule we are implicitly minimizing the MSE of the estimate with respect to the hidden true 
parameter 9. 

Unfortunately, the general solution of equation ^ is intractable, as in principle it requires the 
knowledge of 9, which is the essence of the inference problem. An additional issue with equation (|^ 


^Actually, in a real situation we may not even be able to write such a function due to our imperfect characterization 
or limited knowledge of the detector device. 



















is that, by itself, it does not provide an analytical expression that tells us how to compute 6 in terms 
of /fl Fortunately, there are performance bounds that characterize how far can we be from the 
theoretical solution in equation Q, and even scenarios where the optimal solution can be achieved 
in a closed-form (see Section 3.1 and Appendix [A]) . One of the most significant results in this field 
is the CR minimum variance bound which will be further explained below. 


2.2. The Cramer-Rao bound 

The CR bound offers a performance bound on the variance of the family of unbiased estimators. 
More precisely, let be a collection of independent observations that follow a parametric pmf 

ffj defined on N. The parameters to be estimated from {h,l 2 , will be denoted in general by 

the vector 6 = {9i, 02 ,..., 9m) £ © = Let Tn{I^) : —)• 0 be an unbiased estimatoi]^ of 0, and 

L{I^; 9) = fg{Ii) • fgih) • • • feil-n) be the likelihood of the observation G N" given 9 e 0. Then, 
the CR bound ( Rao|[1945 Cramer||1946 ) establishes that if: 




then, 

Var{TniI'")i) > [l0{n)~^]i,i, 

where 1g{n) is the Fisher information matrix given by: 

jdlnL{I-;9) dlnL{I^-9) 


d9i 


804 


Vi, j E {1,... ,m} . 


( 7 ) 


( 8 ) 


(9) 


In particular, for the scalar case (m = 1), we have that for all 0 E 0: 


min Var(Tn(I^)) >Ie(n) ^ =Ejn^f7i 
r„(-)er" ® 


/dlnL{r-,9)^ ^ 

V M 


( 10 ) 


where is the collection of unbiased estimators and F ^ fjf. 


Returning to our problem in Section 2.1, Mendez et al. (2013 2014) have characterized and 


analyzed the CR bound for the isolated problem of astrometry and photometry, respectively, as 
well as the joint problem of photometry and astrometry. Particularly, we highlight the following 
results, which will be used later on: 


®In many cases one can propose an “educated guess” for 6 (e.g., a function of some sort) which is trained with a 
(“good”) subset of the data itself (thus approximately removing the ambiguity that the true parameter 6 is, in fact, 
unknown). This heuristic approach is adopted, e.g., in the SDSS pipeline trough the use of the KL transform, which 


is a good approximation to a matched filter (also known as the “North filter”, see e.g., Gray & Davisson (2010l). 
®In the sense that, for all 0 G 0, {'rn(l")} = 6- 
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Proposition 1 (Mendez et al. 1(2014 , PP- 800)) Let us assume that Xc G M is fixed and known, 
and we want to estimate F (fixed but unknown) from /"■ ~ f{xcF) equation (5). In this scalar 
parametric context, the Fisher information is given by: 




gfixc 


^ Fgfixc) +Bi 


( 11 ) 


which from equation (10) induces a MV bound for the photometry estimation problem. On the 


other hand, if F € M"*" is fixed and known, and we want to estimate Xc (fixed but unknown) from 
in equation then the Fisher information is given by: 

2 




, dgi(xc) 

dxc 


1 


Fgi{Xc) + Bi ^CR 


( 12 ) 


which from equation (10) induces a MV bound for the astrometric estimation problem, and where 
a‘(jj^ = Ixfin)~^ denotes the (astrometric) CR bound. 


At this point it is relevant to study if there is any practical estimator that achieves the CR 
bound presented in equations © and ( |12[ ) for the photometry and astrometry problem, respec¬ 
tively. For the photometric case, [Perryman et al. (1989, their Appendix A) has shown that the 
classical LS estimator is near-optimal, in the sense that its variance is close to the CR bound for a 
wide range of experimental regimes, and furthermore, in the low SNR regime, when Fgfixc) <C Bi, 


its variance (determined in closed-form) asymptotically achieves the MV bound Fp{n) ^ in equa¬ 
tion ©• This is a formal justification for the goodness of the LS as a method for doing isolated 
photometry in the setting presented in Section [2.1[ An equivalent analysis has not been conducted 
for the astrometric problem, which is the focus of the next section of this work. 


3. Achievability analysis of the Cramer-Rao bound in astrometry 


We first evaluate if the CR bound for the astrometric problem, from equation (12), can 
possibly be achieved by any unbiased estimator. Then, we focus on the widely used LS estimation 
approach, to evaluate its performance in comparison with the astrometric MV bound presented in 
Proposition [TJ 


3.1. Non achievability 

Concerning achievability, we demonstrate that for astrometry (i.e., assuming F is known) 
there is no estimator that achieves the CR bound in any observational regime The log-likelihood 


^^However, in Section]^ we demonstrate that in the low SNR limit the LS estimator can asymptotically approach 
the CR bound. 
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function associated to equation ([^ in this case is given by: 


InL(/”;Xc) = InAi(xc) - Aj(xc) - In/*!), (13) 

i=l 

and we have the following result: 

Proposition 2 For any fixed and unknown parameter Xc G M, and any unbiased estimator r-n 

Var{Tn{P)) > alp., (14) 

where /"■ follows a Poisson pmf (f\[x^) = fx^ hereafter, to shorten notation) from equation 
(The proof is presented in Appendix\^. 

The non-achievability condition imposed by Proposition [^supports the adoption of alternative 
criteria for position estimation, being ML and the classical LS two of the most commonly adopted 
approaches. The ML estimate of the position is obtained through the following rule: 


= argmaxlnL(/”; a) 
aSM 


(15) 


where “arg max” represents the argument that maximizes the expression, while a is a generic vari¬ 
able representing the astrometric position. Imposing the hrst order condition on this optimization 

dXviLil^'x 1 

problem, it reduces to satisfying the condition - } ’ = 0, and, consequently, we can work 


with the general expression given by equation (A2). We note that a well-known statistical result 


indicates that in the case of independent and identically distributed samples the ML approach 
is asymptotically unbiased and efficient (i.e., it achieves the CR bound when the number of ob¬ 
servations goes to infinity (Kay 1993[ Chapter 7.5)). However for the (still independent) but 
non-identically distributed setting of astrometry described by equation (§, this asymptotic result, 
to the best of our knowledge, has not been proven, and remains an open problem. 


On the other hand, a version of the LS estimator (given the model presented in Section 2.1) 
corresponds to the solution of: 


tls( 7”) = argmin^ (Ij - , 

oSM 


i=l 


= J{I",a) 


with Aj(a) = Fgfia) + Bi and where gfi-) is given by equation ([s 


11 


(16) 


^^Note that if we had assumed a detection process subject to a purely Gaussian noise with an rms of a, then 
the probability mass function for each individual observation would have been given by fx{x) = e~ ■ 
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In a previous paper ( Mendez et al.|[^013 Section 5), we have carried out numerical simulations 
using equations (15) and (16), and have demonstrated that both approaches are reasonable. How¬ 
ever, an inspection of Mendez et al. (2013, Table 3) suggests that the LS method exhibits a loss 
of optimality at high-SNR in comparison with either the (Poisson variance-) weighted LS or the 
ML method. This motivates a deeper study of the LS method, to properly understand its behavior 
and limitations in terms of its MSE and possible statistical bias, which is the focus of the following 
Section. 


3.2. Bounding the performance of the Least-Squares estimator 


The solution to equation (16) is non-linear (see Figure]^, and it does not have a closed-form 
expression. Consequently, a number of iterative approaches have been adopted (see, e.g.. Stetson 


(1987); Mighell (2005)) to solve or approximate Hence as is implicit, it is not 


possible to compute its mean, its variance, nor its estimation error directly. We also note that, since 
we will be mainly analyzing the behavior of rL 5 (/”), none of the caveats concerning the properness 
of the likelihood function (equation ([^) raised in Section 2.1 are relevant in what follows, except for 
what concerns the adequacy of equation which we take as a valid description of the underlying 
flux distribution. 

The problem of computing the MSE of an estimator that is the solution of an optimization 


problem has been recently addressed by So et al. (2013) using a general framework. Their basic idea 


was to provide sufficient conditions on the objective function, in our case J(/”', a), to derive a good 
approximation for |('^Ls(7"') — Xc)^|. Based on this idea, we provide below a refined result 

(specialized to our astrometry problem), which relaxes one of the idealized assumptions proposed in 
So et al. (2013, their equation (5)), and which is not strictly satisfied in our problem (see Remark]^ 


in Section 3.3). As a consequence, our result offers upper and lower bounds for the bias and MSE 


of T£S'(I"'), respectively. 


Theorem 1 Let us consider a fixed and unknown parameter x 
addition, let us define the residual random variablfi'^ 
there exists 6 G (0,1) such that ¥{W{I^ ,Xc) G (—5, <5)7^= 1, then: 

- Xc\ < e((5), 


. G M, and that I" ~ fx,,- In 

Win rvi = Jf 

yy(I ,a) - 




(17) 

(18) 


In this case the log-likelihood function (i.e., the equivalent of equation (13l) would be given by lnL(7”;Xc) = 
-nlnjv^o-) - (T - \i{xc)fi. Therefore, in this scenario, finding the maximum of the log-likelihood would 


be the same as finding the minimum of the LS, as in equation (16l. This is a well-established result, described in 
many statistical books. 

^^As a short hand notation: J'{n,Xc) = I 


,,and J"(7",x,) = 
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where 


2 , ^ 

CFjqin) = -^ 

1-5' 


(The proof is presented in Appendix . 


(19) 


( 20 ) 


3.3. Analysis and interpretation of Theorem 


Remark 1 Theorem^is obtained under a bounded condition (with probability one) over the ran¬ 
dom variable W{I^,Xc)- To verify whether this condition is actually met, it is therefore important 


to derive an explicit expression for W{I^,Xc)- Starting from equation (16) it follows that: 


r(r,a) = 2Y. 


i=l _ 


dXi{a) 

da 


+ (Ai(a) — 1%) ■ 


df\i{ 


a} 


d?a 


and, consequently, {J"(/”,Xc)} = ‘2^Ya=i ( ‘^^d!y'^ \a=x„) ■ Therefore: 


W{P,Xc) = Y,iXi{xc)-h)- 


2 = 1 


A"(.t=)/J;(A'(!,))= 

i=i 


( 21 ) 


Then, W{I^,Xc) is not bounded almost surely, since Ii could take any value in N with non-zero 
probability. However, {W{H,Xc)} = 0, and its variance in closed-form is: 


Var{W{H,Xc)) = ^Ai(xc) 


2=1 


(A?(xo))V I 

i=i 


( 22 ) 


From this, we can evaluate how far we are from the hounded assumption of Theorem\^ To do this, 
we can resort to Markov’s inequality {Cover & Thorrias 2006), where ¥{W{H, Xc) i {-p,p)) < 
VarfW{H,Xc))/p^■ Then, for any e G (0,1), we can characterize a critical 5(e) > 0 such that 
F{W{H,Xc) G (—5(e), 5(e))) > 1 —e. Using this result and Theorem^ we can bound the conditional 
bias and conditional MSE of using equations (H and HI, respectively. In Section^ we 

conduct a numerical analysis, where it is shown that the bounded assumption for W{H,Xc) is 
indeed satisfied for a number of important realistic experimental settings in astrometry (with very 
high probability). 
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Remark 2 Concerning the MSE of the LS estimator, equation (18) offers a lower and upper bound 


in terms of a nominal value cr‘j^g{n) (given by equation (19)), and an interval around it. In the 


interesting regime where (5 <C 1 (this regime approaches the ideal case 5 = 0 studied by So et al 


(2013) in which case the variable W{I'^,Xc) becomes deterministic), we have that tls{') is an 
unbiased estimator, as shown by equation 0. and, furthermore: 

Var{TLs{n) = = als{n) > 


from equations (IS) and (10). Thus, it is interesting to provide an explicit expression for a\g{n) 
which will be valid for the MSE of the LS method in this regime. First we note that, J'{I^,Xc) = 
2 • (l]r=i and therefore: 


(x,))2 = 4 • ^ ^ {liff - hXj - IjXi + AiA 


j=l j=l,j^i 
n 


dXi dXj 
^ dxc dxr 


+4 


J;(/2-2/,A, + A2).(^^ 

i=l ^ 


Therefore | Xc)^} = 4 • Yll(=i which implies that: 




Er=i A.(x.) • (A^(xe))^ Er=i iF9i{xc) + B^) • {g[{xfff 


(Er=i(A'(xc))2) 


2^2 




(23) 


In the next section, we provide a numerical analysis to compare the predictions of equation (23) 


with the CR bound computed through equation (12). We also analyze if this nominal value is 


representative of the performance of the LS estimator. 


Remark 3 (Idealized Low SNR regime) Following the ideal scenario where 5 <C 1, we explore the 
weak signal case in which Fgi(xc) 'C Bi considering a constant background across the pixels, i.e.. 


Bi = B for all i. Then adopting equation (23) we have that: 




B 


F^T:=Mi^c)r^ 


(24) 


On the other hand, from equation (12) we have thatZ^ffn) k. /B X)r=i(5i(^c))^- Remarkably 

in this context, the LS estimator is optimal in the sense that it approaches the CR bound asymp¬ 
totically when a weak signal is observec , This results is consistent with the numerical simulations 
Mendez et al. (2013. Table 3)^ 


in 


^^Recall that, according to Sectionjs^ we have demonstrated that the CR can not be exactly reached in astrometry. 

^“^We note that this asymptotic result can be considered the astrometric counterpart of what has been shown in 
photometry by Perryman et al. 119891, where the LS estimator also approaches the CR bound in the low SNR regime. 
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Remark 4 (Idealized High SNR regime) For the high SNR regime, assuming again that 5 <til, we 
consider the case where Fgiixc) S> Bi for all i. In this case: 




2^2 


. (Er=igK^c)^) 

Y17=i9iixc)gi{xc 


1 -1 


and 


1 -1 


^'^i9'i{xc)f / gi{xc 


i=l 


(25) 


Therefore, in this strong signal scenario, there is no match between the variance of the LS estimator 
and the CR bound, and consequently, we have that cr'j^g(n) > To provide more insight on the 

nature of this performance gap, in the next proposition we offer a closed-form expression for this 
mismatch in the high-resolution scenario where the source is oversampled, and the size of the pixel 
is a small fraction of the width parameter a of the PSF in equation 0. 

Proposition 3 Assuming the idealized high SNR regime, if we have a Gaussian-like PSF and 
Ax fa <C I, then: 


a 


LS 


n 


8 


a. 


CR 


3\/3 


> 1 . 


(26) 


(The proof is presented in Appendix . 


Equation (26) shows that there is a very significant performance gap between the CR bound and 
the MSE of the LS estimator in the high SNR regime. This result should motivate the exploration 
of alternative estimators that could approach more closely the CR bound in this regime. 


4. Numerical analysis of the LS estimator 


In this section we explore the implications of Theorem in astrometry through the use of 
simulated observations. Eirst, we analyze if the bounded condition over W{H,a) adopted in 
Theorem is a valid assumption for the type of settings considered in astronomical observations. 
After that, we compare how efficient is the LS estimator proposed in equation ( |16| ) as a function 
of the SNR and pixel resolution, adopting for that purpose Theorem and Proposition 


To perform our simulations, we adopt some realistic design variables and astronomical observ¬ 
ing conditions to model the problem ( Mendez et al.|20l3 2014). For the PSF, various analytical and 
semi-empirical forms have been proposed, see for instance the ground-based model in King (1971) 
and the space-based model in Bendinelli et al. (1987). For our analysis, we adopt in equation ([s]) 

-I _ ^ ' 

the Gaussian PSF where 4>{x,a) = ^ and where a is the width of the PSF, assumed 

to be known. This PSF has been found to be a good representation for typical astrometric-quality 
ground-based data (Mendez et al. 2010). In terms of nomenclature, the FWHM = 2\/2 In 2 a, 
measured in arcsec, denotes the Full-Width at Half-Maximum (FWHM) parameter, which is an 


overall indicator of the image quality at the observing site (Chromey 


2010 


15 


our simulations we fix the value of FWHM to a constant value. 


However, in many cases the image quality 






























- 15 - 


The background profile, represented by |-Bi, i = 1, n|, is a function of several variables, like 

the wavelength of the observations, the moon phase (which contributes significantly to the diffuse 
sky background), the quality of the observing site, and the specifications of the instrument itself. 
We will consider a uniform background across pixels underneath the PSF, i.e., Bi = B for all i. To 
characterize the magnitude of B, it is important to first mention that the detector does not measure 
photon counts (or, actually, photo-e“) directly, but a discrete variable in Analog to Digital Units 
(ADUs)” of the instrument, which is a linear proportion of the photon counts (Gilliland 1992). 


This linear proportion is characterized by the gain of the instrument G in units of e“/ADU. G is 
just a scaling factor, where we can dehne F = F/G and B = B/G as the brightness of the object 
and background, respectively, in the specihc ADUs of the instrument. Then, the background (in 
ADUs) depends on the pixel size Ax as follows (Mendez et al.||2013): 


B = fs 


^ D + RON^ 

Ax H- — - ADU, 

Ct 


(27) 


where fs is the (diffuse) sky background in ADU/arcsec (if Ax is measured in arcsec), while D and 
RON, both measured in e~, model the dark-current and read-out-noise of the detector on each 


pixel, respectively. Note that the first component in equation (27) is attributed to the site, and 


its effect is proportional to the pixel size. On the other hand, the second component is attributed 
to errors of the integrating device (detector), and it is pixel-size independent. This distinction 
is central when analyzing the performance as a function of the pixel resolution of the array (see 


details in Mendez et al. (2013, Sec. 4)). More important is the fact that in typical ground-based 
astronomical observation, long exposure times are considered, which implies that the background 
is dominated by diffuse light coming from the sky (the first term in the RHS of equation (|27[)), and 


not from the detector (Mendez et al. 2013, Sec. 4). 


For the experimental conditions, we consider the scenario of a ground-based station located at 
a good site with clear atmospheric conditions and the specifications of current science-grade CCDs, 
where fs = 1502.5 ADU/arcsec, D = 0 e~, RON = 5 e~, FWHM = 1 arcsec and G = 2 e“/ADU 
(with these values B = 313 ADU for Ax = 0.2 arcsec using equation (|27[)). In terms of scenarios 
of analysis, we explore different pixel resolutions for the CCD array Ax G [0.1, 0.7] measured in 
arcsec. Note that a change in Ax will be reflected upon the limits of the integral to compute the 
pixel response function gi{xc) (see equation as well as in the calculation of the background level 
per pixel Bi, according to equation (27). Therefore, a change in Ax is not only a design feature of 


the detector device, but it implies also a change in the distribution of the background underneath 
the PSF. The impact of this covariant device-and-atmosphere change in the CR bound is explained 


in detail in Mendez et al. (2013, Section 4, see also their Figure 2). 


changes as a function of position in the field-of-view due to optical distortions, specially important in large focal 
plane arrays. For example, in the case of the SDSS (which consists of 30 CCDs of 2048^ pix^ each, covering 2°.3 on 
the sky at a resolution of 0.396 arcsec/pix), the FWHM may vary up to 15% from center to corner of one detector 
(Lupton et al. (2001 Section 4.1)). The impact on the CR bound of these changes has been discussed in some detail 


by 


Mendez et al. 


(2014 Section 3.4). 
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In our simulations, we also consider different signal strengths F G {1 080, 3 224,20 004, 60 160}[^ 
measured in photo-e“, corresponding to SNR G~ {12,33,120,230} respectivelj|^ Note that in¬ 
creasing F implies increasing the SNR of the problem, which can be approximately measured by 
the ratio F/B = F/Rg On a given detector plus telescope setting, these different SNR scenarios 
can be obtained by changing appropriately the exposure time (open shutter) that generates the 
image. 


4.1. Analyzing the bounded condition over 


To validate how realistic is the bounded assumption over W{I^,q) in our problem, we hrst 
evaluate the variance of W{I'^,Xc) from equation (22), this is presented in Figure for different 
SNR regimes and pixel resolutions in the array. Overall, the magnitudes are very small considering 
the admissible range (0,1) for VF(/”,Xc) stipulated in TheoremAlso, given that VF(/",Xc) has 
zero mean, the bounded condition will happen with high probability. Complementing this, Figurej^ 
presents the critical 5 across different pixel resolutions and SNR regime^] For this, we fix a small 
value of e (= 10“^ in this case), and calculate 5 such that W{I"',Xc) G (—with probability 
1 — e. From the curves obtained, we can say that the bounded assumption is holding (values of 6 
in (0,1)) for a wide range of representative experimental conditions and, consequently, we can use 
Theorem to provide a range on the performance of the LS estimator. Note that the idealized 
condition of J ~ 0 is realized only for the very high SNR regime (strong signals). 


4.2. Performance Analysis of the LS estimator 


We adopt equation (18) which provides an admissible range for the MSE performance of the 
LS estimator. For that we use the critical 5 obtained in Figure These curves for the different 
SNR regimes and pixel resolutions are shown in Figure Following the trend reported in Figure 
the nominal value is a precise indicator for the LS estimator performance for strong signals 
(matching the idealized conditions stated in Remark]^, while on the other hand. Theoremdoes 
not indicate whether is accurate or not for low SNR, as we deviate from the idealized case 
elaborated in Remark Nevertheless, we will see, based on some complementary empirical results 


“These are the same values explored in 


Mendez et al. 


(2013 


Table 3). 


^For a given F and fs there is a weak dependency of SNR on the pixel size Ax, see equation (28) in 


Mendez et al. 


(20131. 

We note that while the ratio F/B can be used as a proxy for SNR, in what follows we have used the exact 
expression to compute this quantity, as given by equation (28) in 


Mendez et al. 


(20131. 


^®These values were computed empirically (from frequency counts) using 5 000 realizations of the random variable 
1T(7”, Xc) for the different SNR regimes and pixel resolutions. 
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reported in what follows, that even for low SNR, the nominal predicts the performance of the 
LS estimator quite well. 


Assuming for a moment the idealized case in which <5 <C 1, we can reduce the performance 
analysis to measuring the gap between the nominal value predicted by Theorem(equation (|19[)), 


and the CR bound in Proposition ll Figure 




shows the relative difference given by e% = 




j2 

^CR 


100. From the figure we can clearly see that, in the low SNR regime, the relative performance 
differences tends to zero and, consequently, the LS estimator approaches the CR bound, and it is 
therefore an efficient estimator. This matches what has been stated in Remark[3l On the other hand 
for high SNR, we observe a performance gap that is non negligible (up to ~ 27% relative difference 
above the CR for F = 60 160 e~, and ss 15% above the CR for F = 20 004 e~ for Ax = 0.2 arcsec). 
This is consistent with what has been argued in Remark]^ Note that in this regime, the idealized 
scenario in which 5 <C 1 is valid (see Figurej^ and, thus, |(ri 5 (/’^) — Xc)^| ~ which 

is not strictly the case for the low SNR regime (although see Figure and the discussion that 
follows). 


To refine the relative performance analysis presented in Figure]^ Figureshows the feasible 
range (predicted by Theoremof performance gap considering the critical 6 obtained in Figure]^ 
We report four cases, from very low to very high SNR regimes, to illustrate the trends. From this 
figure, we can see that the deviations from the nominal value are quite significant for the low SNR 
regime, and that, from this perspective, the range obtained from Theorem[^is not sufficiently small 
to conclude about the goodness of the LS estimator in this context. On the other hand, in the high 
SNR regime, the nominal comparison can be considered quite precise. 


The results of the previous paragraph motivate an empirical analysis to estimate the per¬ 
formance of the LS estimator empirically from the data, with the goal of resolving the low SNR 
regime illustrated in Figure For this purpose, 1000 realizations were considered for all the SNR 
regimes and pixel sizes, and the performance of the LS estimator was computed using the empir¬ 
ical MSE. We used a large number of samples to guarantee convergence to the true MSE error as 
a consequence of the law of large numbers (Gray & Davisson 2010). Remarkably, we observe in 
all cases that the estimated performance matches quite tightly the nominal characterized by 
Theorem We illustrate this in Eigure|^ which considers the most critical low SNR regime 
Consequently, from this numerical analysis, we can resolve the ambiguity present in the low SNR 
regime, and conclude that the comparison with the nominal result reported in Figure and the 
derived conclusion about the LS estimator in the low and high SNR regimes, can be considered 
valid. Our simulations also show that the LS estimator is unbiased. Overall, these results suggests 
that Theorem could be improved, perhaps by imposing milder sufficient conditions, in order to 
prove that a\g is indeed a precise indicator of the MSE of the LS estimator at any SNR regime. 

Eor completeness, we show in Eigure[^the behavior of the log-likelihood function (computed 


^^Similar results were obtained in the other regimes, not reported here for the sake of brevity. 
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using equation ([I^) and the LS function J{I'^,Xc) (computed using equation dlh])) for our two 
extreme SNR cases of our numerical simulations, namely SNR= 12 and SNR= 230. In these figures 
the true astrometric position is at Xc = 80 arcsec (= 400 pix with Ax = 0.2 arcsec). These figures 
clearly show (particularly the one at low SNR) the non-linear nature of both objective functions. 


5. Summary and Final Remarks 


Our work provides results to characterize the performance of the widely used LS estimator as 
applied to the problem of astrometry when derived from digital discrete array detectors. The main 
result (Theoremj^ provides in closed-form a nominal value (c’'|g(n)), and a range around it, for the 
MSE of the LS estimator as well as its bias. From the predicted nominal value, we analyzed how 
efficient is the LS estimator in comparison with the MV CR bound. In particular, we show that the 
LS estimator is efficient in the regime of low SNR (a point source with a weak signal), in the sense 
that it approximates very closely the CR bound. On the other hand, we show that at high SNR 
there is a significant gap in the performance of the LS estimator with respect to the MV bound. We 
believe that this sub-optimal behavior is caused by the Poissonian nature of the detection process, 
in which the variance per pixel increases as the signal itself. Since the LS method is very sensitive 
to outliers, the large excursions caused by the large pixel intensity variance at high SNR make the 
LS method less efficient (from the point of view of its MSE), than allowed by the CR bound. These 
performance analyses complement and match what has been observed in photometric estimation, 
where only in the low SNR regime the LS estimator has been shown to asymptotically achieve the 
CR bound. 


While our results are valid for an idealized linear (one-dimensional) array detector where intra¬ 
pixel response changes are neglected, and where flat-fielding is achieved with very high accuracy, our 
findings should motivate the exploration of alternative estimators in the high SNR observational 
regime. Regarding this last point, we note that an inspection of Mendez et al. (2013, Table 3) 
suggests that either a (Poisson variance-) weighted LS or a ML approach do not exhibit this loss 
of optimality at high-SNR, and should be preferred to the unweighted LS analyzed in this paper. 
This effect is clearly illustrated in Figure where we present a comparison between the standard 
deviation of the LS and ML methods derived from numerical simulations in the very high SNR 
regime (where the gap between CR and the LS is most significant). Motivated by these results, a 
detailed analysis of the ML method will be presented in a forthcoming paper. 
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A. Proof of Proposition 

We use the well-known fact that the CR bound is achieved by an unbiased estimator, if and 
only if, the following decomposition holds (Kendall et al.||1999 pp. 12) 


dlnL{r-9) 


= A{e,n)-{f{n-e), 


(Al) 


where 9) is the likelihood of the observation I'^ G N"’ given 0, A(0, n) is a function of 9 and n 
alone (i.e., it does dot depend on the data), while /(J*^) is a function of the data exclusively (i.e., 
it does not depend on the parameter). Furthermore, if the achievability condition in equation (Al) 
is satisfied, then A{9,n) = and /(/’*) is an unbiased estimator of 9 that achieves the CR 

bound. 


The proof follows by contradiction, assuming that equation (Al) holds. First using equa¬ 
tion (j^, we have that: 


din L(I’^;Xc 
dxr 




d\i{xc) d\i{xc 


dxc 


dxr 


E 


Aj(xc) 


dXjjxc) 

dXr 


(A2) 


the last equality comes from the fact that Y17=i = 1 from the assumption in equation (j^. 

Then replacing equations (A2) and (12) in equation (Al): 


fin = E 


. dgijxc) 

F9i(xc)+Bi dXc 




1 


Fgi{xc)+Bi 


dgijxc) 
dxc 


+ Xc 


(A3) 


which contradicts the assumption that f{I^) should be a function of the data alone. Furthermore, 
if we consider the extreme high SNR regime, where Fgi{xc) ^ Bi for all i, and the low SNR regime, 
where Fgi{xc) <C Bi for all i, it follows that: 


i=l 


dhigi{xc) 

dxc 


and 




li dgi{xc) 


1 Bi 
1=1 ‘ 


dXr 


^E 

i=l 


^E 

i=l 


dgijxc 

dxc 


dgijxc) 

dXr 


/giiXc 


/Bi 


-1 


+ Xc 


+ Xc 


(A4) 


(A5) 


respectively. Therefore a contradiction remains even in these extreme SNR regimes. 
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B. Proof of Theorem [T] 


The approach of So et al. (2013) uses the fact that the objective function J{I"’,a) in equa¬ 


tion (16) is two times differentiable, which is satisfied in our context. As a short-hand, if we denote 

l ocal o ptimum 
(2013) is that 


by Xc the LS estimator solution then the first order necessary condi tion for a l ocal o ptimum 

requires that J'{I'^,Xc) = | a=xc — 0- The other key assumption in 


So et al. 


Xc is in a close neighborhood of the true value Xc- In our case, this has to do with the quality of 
the pixel-based data used for the inference, which we assume it offers a good estimation of the po¬ 


sition (see, e.g.. King (1983); Stone (1989)). Then using a first order Taylor expansion of J'(/"',Xc) 

I): 

(Bl) 


around Xc, the following key approximation can be adopted (So et al. 2013, their equation (4)): 


0 = J'{P, Xc) = Xc) + (Xc - Xc) ■ Xc) ^ {Xc - Xc) = 


t// jn 


// / jn 


j"{i^,xcy 


where J"{I'^,Xc) = \a=xc- If we consider I" ~ then from equation (Bl): 

J'{r,xc) 


Xc - TLsil"") = 


J”{I^,Xc)' 


(B2) 


The second step in the approximation proposed by |So et al, 


(2013) is to bound by 




{j'qj" X )} ■ Tor that we introduce the residual variable W a) where Xc) = Xc)} (14 


W (/”■, Xc)). Using the fact that W (I”, a) is bounded almost surely (see Remark]^ and Section 4.1): 


j'(/-,Xc) J'(/^Xc) 


J'(/”,Xc) 

r 1 1 

E7"~/.,{J"(/",Xc)} J"(/-,Xc) 


{T"(-I’",Xc)} 

1 4 IU(/’',Xc)J 


< 


J'{I^,Xc) 


< 


^i^^uAJ"ii^,xc)} 

\J'{I",Xc)\ 


max 

w€.{—S,5) 


1 - 


1 


1 + w 


T/"~A.{T"(/",Xc)} 1-<5’ 


(B3) 


the last step uses the fact that {J"(/’^, Xc)} > 0 (see Remark [y. On the other hand, 

Jensen’s inequality ([Cover fc: Thomas 2006) guarantees that: 


T/"~A.{T'(/",Xc)} 






[ J'jl'^.Xc) 

U"(/”,Xc) 




J'{r,xc) 


J'{I^,Xc) 


< 


T/"~A.{T"(/",Xc)} J"(/-,Xc) 

T/-~A.{|T'(/”,Xc)|} <5 




(B4) 


where the last inequality comes from equation (B3). Then we use that J'(I”, Xc) = 2 -X)r=i (AlI^c) ~ 
li) ■ and consequently {J'(/",Xc)} = 0. Then from equations (B4) and (B2), we 

have that: 

I TIT r ^ {|T (.!"')2:c)|} 5 

|xc - {tls{I )}| < - r in nr,. ^ - 1 ’ (^ 5 ) 


{«/"(/", Xc)} 1-<I’ 


which leads to equation (0. 
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Concerning the MSE, from the hypothesis on W{I'^,Xc) we have that: 




< 




Xc)f (1 + 5)2 Xe) ) Xc)f (1 - 5)2 ’ 


< 




(B6) 


almost surely. Then taking the expected value in equation (B6) and using equation (B2) for the 
central term, it follows that: 


(1 + 5)2 

which concludes the result. 


< 




< 


(1-5)2’ 

(B7) 
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C. Proof of Proposition 


^. I — 

Recalling from equation (3) that gi{xc) = f ‘ aI ~ Xc,cr)dx and assuming a Gaussian 

Xi ^ 


PSF of the form (p{x,a) = exp ^ 2 ^) (see Section 

hypothesis of small pixel (Ax/u <C 1), it is possible to state that: 

gi{xc) ~ 4>ixi - Xc, cr) ■ Ax, 

and then we have that: 


4) by the mean value theorem and the 

(Cl) 


dgi(xc) (xi-Xc) 1 


dXr 


o 


exp 


-{Xi - Xcf 
2a^ 


• Ax 


{Xi - Xc) ,, . . 

- 5 - (plxi - Xc,cr) ■ Ax. 


(7^ 


(C2) 


With the above approximation we have that: 


E 

2=1 


dgijxc) 

dXr 


giiXc 


E 


{Xi - Xc) 


2 = 1 
Ax^ 




(j){xi — Xc, cr) ■ Ax ) • (j){xi — Xc, cr) ■ Ax 


LAX \2 X 

—^-• > (Xi — Xc) , —exp 


-{Xi - Xcf 


\ 


• Ax. (C3) 


The term inside the summation in equation (C3) can be approximated by an integral due to the 


small-pixel hypothesis. Assuming that the source is well sampled by the detector (see Section 2.1 
equation Q) we can obtain that: 

/ 


E 

2=1 


dgijxc) 

dxc 


giixc) 


Ax^ 


2\/37r(T6 7_oo 


-1-00 1 

(x - Xc)^ 




exp 


V3 


-(x - Xc)" 


Ax^ cr^ Ax^ 

2\/37r(T® 3 6\/37riT^ ’ 


dx (C4) 


(C5) 


where equation (C5) follows from the fact that the term inside the integral in equation (C4) corre- 
' 2 ' 
sponds to the second moment of a normal random variable of mean Xc and variance 


By the same set of arguments used to approximate Y7i=i gi{xc) in equation (C5), we 

have that: 


E 

2=1 


dgijxc) 

dXr 


E 

2=1 

Ax 


(Xi - Xc) 




(f){xi — Xc, cr) ■ Ax 






-(Xj - Xc)^ 


V 


• Ax 
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Ax 


Ax 


/ + 00 1 


/ 


exp 


72 


-{X - Xc)" 


V 


dx 


a 


Ax 


2 y^(T^ 


4y^cr^ 


(C 6 ) 


(C7) 


where equation (C7) follows from the fact that the term inside the integral in equation (C 6 ) is the 
' 2 ' 
second moment of a normal random variable of mean Xc and variance 


Finally, for 


n ( dgi{xc) 


dxc 


9i{^c) 


we proceed again in the same way, namely: 


E 

1 = 1 


dgijxc) 

dXr 


gi{Xc 


E 

2 = 1 


{Xi - Xc) 




(j)(xi — Xc, cr) ■ Ax 
1 


4>{xi — Xc, cr) ■ Ax 


a 


2 = 1 ^ 


-{xi - Xc)^ 


Ax 


a 


1 f+^{x-xc)^ 

-4 / — 


' —00 


\/^i 


TTCr 


-{x - Xcf 
2^2 


dx 


1 ^ 2^1 
a 2 ’ 


(C 8 ) 

(C9) 


where equation (C9) follows from the fact that the term inside the integral in equation (C 8 ) corre¬ 
sponds to the second moment of a normal random variable of mean Xc and variance cr 2 . 


Then adopting equations (C5), (C7) and (C9) in equations (23) and ( 12 ), respectively, and 


assuming a uniform background underneath the PSF, we have that: 
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In the last step in equation (CIO) and the last step in equation (Cll), we have used the assumption 
that F ^ B. We note that the last step in equation (Cll) corresponds to the second line of 
equation (45) in Mendez et al. (2013). 
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Fig. 1.— Variance of the residual random variable VF(/"',Xc) (dimensionless) as given by equa¬ 
tion (22), as a function of the pixel resolution Ax arcsec for different realistic SNR scenarios 
(function of F) encountered in ground-based astronomical observations. Since the admissible range 
for 5 is the interval (0,1), the small computed values indicates that the bounded assumption in 
Theorem 1 can be considered as valid under these conditions. 
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Fig. 2.— Numerical computation of the critical 6 (dimensionless) such that P{W{I^,Xc) G 
(—(5,5)) > 1 — 10“^. In all the scenarios (SNR, and Ax), 5 000 realizations of the random variable 
W{I^,Xc) are used to estimate the probability distribution for W{I^,Xc), and 5, from frequency 
counts. As S decreases, we have a smaller bias (see equations 0 and (|20[)) and a narrower range 
for the MSE of the LS estimator (equation ([I^). 
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Fig. 3.— Range of the MSE performance (in milli-arcsec^ = mas^) for the LS method in astrometry 
predicted by Theorem]^ (equations and ([T^). Results are reported for different representative 
values of F and across different pixel sizes: (From top-left to bottom-right) F = 1080 e“; F = 
3 224 e"; F = 20 004 e"; F = 60160 e". 
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Fig. 4.— Relative performance differences between in Theorem (equation (|l9[)) and the CR 
bound in Proposition (equation (12)). Results are reported for different SNR and pixel 
sizes. A significant performance gap between the LS technique and the CR bound is found for 
FWHM/Ax < 1 (good sampling of the PSF) at high SNR, indicating that, in this regime, the LS 
method is sub-optimal, in agreement with Proposition]^ (see also equation (26)). This gap becomes 
monotonically smaller as the SNR decreases. 
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Fig. 

stipulated in Theorem for the LS estimator (equations (18) and (19)) and the CR bound a\^ 
in Proposition]^ (equation (12)). Results are reported for different F and across different pixel 
sizes: (From top-left to bottom-right) F = 1080 e~] F = 3 224 e~] F = 20 004 e~; F = 60 160 e“. 
The 0% level corresponds to having achieved the CR bound. Note that as the SNR decreases, 
the bias (see equations © and ( p0| ), and the possible range for the MSE of the LS method (see 
equations ([T^) increase (see also Figure]^. 
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Fig. 6.— Comparison between the nominal value the performance range ^) 


stipulated in Theorem|^ and the empirical estimation ofEjn^f^^{(TLs(-^’^) — Xc) } from simulations 
for a low SNR regime of F = 1 080 e“. The fact that the simulations follow closely the nominal 


value, even at low SNR, justifies the use of a^s as given by equations (19) and (23) as a benchmark 
of the LS method at any SNR. 
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Fig. 7.— Example behavior of the likelihood function lnL(/”,a:c) (upper row) and LS function 
J(/", Xc) (lower row) for two particular realizations with low SNR= 12 (left column) and high SNR= 
230 (right column), for a target centered at Xc = 80 arcsec (for a pixel size of Ax = 0.2 arcsec). All 
the other parameters are those described in Section]^ The figures show (particularly the ones at 
low SNR) the non-linear nature of both objective functions. 


























35 - 



Fig. 8.— Comparison of the y/MSE derived from numerical simulations using the ML (open circles, 
equation p^) and the LS method (open squares, equation (16)) for a high SNR= 230 (see, e.g., 
right column of Figure]^, where the optimality loss (performance gap, ~15% in this case) of the 
LS method (Proposition in this regime is clearly seen. The solid line is the nominal (TLsin) 
derived from our theorem (equation (19)), while the dashed line is the CR limit, (Tcr, given by 
equation (12). As we have shown (Section ( |3.1[ )), the CR limit can not be reached in our astrometric 
setting, but our ML simulations (open circles) show that they can follow very closely this limit (see 
also Mendez et al. (2013, Table 3)). A detailed analytical study of the optimality of the ML method 
will be presented in a forthcoming paper. 







































































































